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Abstract 

Methodological development for the inference of gene regulatory networks from transcriptomic data is an 
active and important research area. Several approaches have been proposed to infer relationships among genes 
from observational steady-state expression data alone, mainly based on the use of graphical Gaussian models. 
However, these methods rely on the estimation of partial correlations and are only able to provide undirected 
graphs that cannot highlight causal relationships among genes. A major upcoming challenge is to jointly analyze 
observational transcriptomic data and intervention data obtained by performing knock-out or knock-down exper- 
iments in order to uncover causal gene regulatory relationships. To this end, in this technical note we present 
an explicit formula for the likelihood function for any complex intervention design in the context of Gaussian 
Bayesian networks, as well as its analytical maximization. This allows a direct calculation of the causal effects 
for known graph structure. We also show how to obtain the Fisher information in this context, which will be 
extremely useful for the choice of optimal intervention designs in the future. 
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1 Introduction 

Inference of gene regulatory networks from transcriptomic data has been a wide research area in recent years. 
Several approaches have been proposed, mainly based on the use of graphical Gaussian models [ 1 1 . These methods, 
however, rely on the estimation of partial correlations and provide undirected graphs that cannot highlight the 
causal relationships among genes. Biihlmann et al. (TJIU recently proposed a method to predict causal effects from 
observational data alone in the context of Gaussian Bayesian networks (GBN). In this method, the PC algorithm 
J3] is first applied to find the associated complete partially directed acyclic graph (CPDAG) among the graphs 
belonging to the corresponding equivalence class. Then, intervention calculus [5] is performed to estimate bounds 
for total causal effects based on each directed acyclic graph (DAG) in the equivalence class. 

If knock-out or knock-down experiments are available, however, it is valuable to perform causal network infer- 
ence from a mixture of observational and intervention data. One approach has been proposed to do so [6 1, based on 
a simple comparison of observed gene expression values to the expression under intervention; the underlying idea 
is that if gene Y is regulated by gene X, then its expression value under a knock-out of gene X will be different 
from the value in a wild type experiment. We note that this method provided the best network estimation in the 
DREAM4 challenge, and has the advantage of being very fast to compute without imposing a restriction on the 
acyclicity of the graph. It does, however, require an intervention experiment to be performed for each gene, which 
can be unrealistic for real applications given the cost and time typically involved for knock-out experiments. In 
addition, although it is well-suited to the inference of the structure of the graph, it tends to be imprecise for the 
estimation of the strength of the interactions between genes. 

The aim of this technical note is to propose an explicit calculation of the likelihood function for complex 
intervention designs, including both observational and intervention data, in the context of GBNs. This calculation 
makes use of the full set of available information available, does not require an intervention for each gene, and is 
able to deal with multiple interventions (e.g., a double gene knock-out experiment). For an known graph structure, 
we present here the likelihood calculation for observational data only, as well as for any intervention design. We 
also provide the analytical first order derivatives which allow a direct estimation of the graph structure as well as 
the causal effects. Finally, we give the Fisher information, which is not trivial to derive, and will be extremely 
useful in the future for the choice of optimal intervention designs. 
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The rest of this technical note is organized as follows. In Section 2, we define the model and set up a toy 
example for illustrative purposes. In Sections 3 and 4, we define the likelihood function, maximum likelihood 
estimators, and Fisher information in the case where only observational data are available and in the case where 
a mixture of intervention and observational data are available, respectively. In Section 5, we provide a brief 
discussion and conclusion. 



2 Model definition 

2.1 Definition 

We consider the set Xx, 1 = {1, . . . ,p} a set of p Gaussian random variables defined by: 

Xj=mj+ ^2 Wi,jXi+£j with e 4 ~ j\f(0, cx 2 ). (1) 

iepa(j) 

We assume that the ej are independent, and that i E pa(j') => i < j; this assumption is equivalent to assuming that 
the directed graph obtained using the parental relationships is acyclic. Given the parental structure of the graph, the 
model parameters are 9 = (to, a, w) where Wij is nonzero only on the edge set (i, j) 6 £ = {i e pa(j'), j £ X}. 
It is easy to see that this model is equivalent to Xx ~ A/"(/x; £), with: 

fi = to,L and S = L T diag(cr 2 )L = cr|L T eJejL 

j 

where ej is a null row- vector except for the its j th term which is equal to 1, and where L = (I — W) _1 = 
I + W + . . . + W p_1 with W = {wi t j) it j e i. Note that the nilpotence of W is due to the fact that w i;j = for 
all i ^ j. 

2.2 A toy example 

We consider the particular case where p = 3, pa(l) = 0, pa(2) = {1}, pa(3) = {1, 2}; the true values of the 
parameters are set to be m* = (0.5 1.2 0.7); a* = (0.3 1.1 0.6); w^ 2 = -0.8, w* 13 = 0.9, and w^ 3 = 0.5. We 
hence have 

-0.8 0.9 \ / 1.0 -0.8 0.5 

W=| 0.5 and L = I + W + W 2 = 1.0 0.5 
/ \ 1.0 

and observed data can be generated through Xi-s ~ A/"(/i; S) with: 

0.090 -0.0720 0.045 
H= (0.5 0.8 1.55) and S=( -0.072 1.2676 0.569 



0.045 0.5690 0.685 



3 Observational data 



3.1 Likelihood 



The log-likelihood of the model described in Equation (Q}, given N observations x k = (x\, . . . , xz) (1 ^ k ^ N), 
may be written as follows: 

e(m, a, w) = log(27r) - N £ logfo) - \ ^ \ YM ~ ~ (2) 

1 3 j °i k 

Proof. For all fc, let us define A k = (x k - mL)S _1 (x fc - toL) t . Since IT 1 = (I - W)diag(l/a 2 )(I - W) T 
we obtain: 

A k = ^i^(x fe (I- W) -m)eje 3 (x k (l- W) - m) T 
3 a i 
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We now consider the derivative of i with respect to m: 

— — [m, a, w) = — 7; > (X; — x We, — m, . 

am,- 07 ' J 

The maximization of ^(m, cr, w) in m for a fixed u> hence does not depend on a and is given by: 

1 N 

fe=l 

By replacing with this formula in Equation (O we get an expression of the likelihood free of the parameter m: 

%, w) = log(27r) - iV £ Iogfa) - ~ ^ 4 - y fc WeJ) 2 (3) 

where for all k,j we have: 

3.2 Maximum likelihood estimator 
3.2.1 Derivatives with respect to w 

The derivatives of I defined in Equation (0 with respect to w are as follows: 
Proof. We first note the following: 

J 1/ fc 1 ■/ ■ J 

y ! 3 

As such, the maximization of i?(cr, w) in u> can be done independently from cr by solving for all (i, j) 6 £: 

iV AT 



E^ feWe J = E^- 



fc=i fc=i 

Hence using 

W = E m'j'e^ej, => y k Wej = E JW* 
{i',j')d£ i',(i',j)e£ 
we find that i£> is solution of the following linear system: 

N N 

E * E Vi Vi' = E Vi fOT a11 (*' 3 ) 6 £ • 

i',(i',j)e£ fe=l fe=l 



3.2.2 Derivatives with respect to cr 

The derivatives of I defined in Equation (0 with respect to cr are: 
The maximization of t(cr, w) in cr when m is fixed is thus given by: 
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3.2.3 Toy example (continued) 

Here is a sample of size N — 5 drawn from our toy example model: 

/ 1.1025540 -0.2652622 1.957083 \ 

0.6721755 0.4286717 1.605024 

0.3455340 2.8835932 1.932982 

0.4139627 1.0847936 1.250889 

\ 0.2844364 1.0490652 1.446954 ) 



where after centering we get: 



/ 0.5388215 -1.30143445 0.31849676 \ 

0.1084430 -0.60750058 -0.03356279 

-0.2181985 1.84742088 0.29439582 

-0.1497698 0.04862127 -0.38769719 

\ -0.2792962 0.01289289 -0.19163260 / 



We obtain w by solving: 



0.4501364 0.0000000 0.000000 \ / w h2 
0.0000000 0.4501364 -1.181107 wi' 3 
0.0000000 -1.1811075 5.478283 / V w 23 



-1.1811075 
0.2153241 
0.1284387 



which gives (reference value in parentheses): 

u>i, 2 = -2.6238878 (w*^ 2 = -0.8) w h3 = 1.2430964 « 3 = 0.9) w 2 , 3 = 0.2914543 (w^ 3 = 0.5). 
Finally, we then obtain: 

<7i = 0.3000455 (a* x = 0.3) a 2 = 0.6898100 (a 2 = 1.1) <r 3 = 0.1193022 (a 3 = 0.6) 
and the nuisance parameter: 

mi = 0.5637325 (ml = 0.5) rh 2 = 2.5153432 (m 2 = 0.2) m 3 = 0.6358156 (ml = 0.7). 



3.3 Fisher information 
3.3.1 Hessian of £ 

The (non-zero) second order derivatives of I are given by 

dH 



dwijdwi> j 



(a, w) 



3 k 



da 



2 (<7,W) 



(o-, «o = - 4 Yl - y kWe 

3 k 



dwijdijj 



J -Ik 



3.3.2 Distribution of y k 

We can rewrite Equation (|4) as: 



N - 1 , 

x k - 

N N 



which is obviously a Gaussian vector with expectation: 

k'^k 



N 



k'^k 

N -1 iV-1 



and variance: 



TV TV 
AT- 1 7V-1 



/x = 



iV 2 



N 2 



-S = 



iV 



It is therefore easy to establish that: 



and y fc - y k W ~JV 0; 



AT- 1 



A^ 



diag(cr 2 ) 



(5) 



Note that as a consequence of this, it is easy to prove that E[<t 2 ] = (AT — 1) /Naj which means that this estimator 
is (slightly) biased. 
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3.3.3 Information 

The Fisher information matrix I(cr, w) can therefore be written as: 

d 2 l 



e Wij l{a,w)e a . = -E 



e m J(a, w)e Wil j 

d 2 e 



-E 



(a, w) 



(a, w) 



9 "^M' 



3 k 



2(N- 1) 



45>E [(/)V] (I-W)eJ = ^l_^Z eiS(I _ w)e J = 



3 k 



because eiS(I — W)eJ = <r 2 eiL T eJ = <j 2 ejLef and L is upper triangular. And finally: 



e aj I(a,w)el. 



-E 



d 2 l 



dwijduj 



(a,w) 



^E[(^-/WeJ) 2 



N _ 3(N- 1) AT _ 2A-3 



3.3.4 Toy example (continued) 

We present here the inverse Fisher information matrix in the particular case of our toy example model. Due to the 
block-wise nature of I(<r, w), the Cramer-Rao lower bound on the covariance matrix is given by blocks: 

2 

(N - l)Var(wi 2 ) = = 13 - 444 
Si i 



( Si,i 


Si, 2 




/ 4.1904132 


0.2380165 


V S 2,l 


^2,2 




■ 0.2380165 


0.2975207 



(A - l)Cov(l&i ) 3,t&2,3) = o\ 

{2N - 3)Cov(<7i,(72, <r 3 ) = diag(cr 2 ) = diag(0.09, 1.21, 0.36). 

In the particular case where N = 200, the standard deviations corresponding to the Cramer-Rao bounds for 

6 = (wi.2, wi,3,wi,3,<7i,<7 2 ,<7 3 ) are: 

sd CR (<9) = ( 0.25992311 0.14511152 0.03866625 0.01505657 0.05520742 0.03011314 ) 
while the empirical standard-deviation (sample size 2000) are: 

sd emp (<9) = ( 0.26110572 0.14847838 0.03873625 0.01476351 0.05424152 0.02952547 ) . 
The empirical mean is: 

mean emp (l9) = ( -0.8024838 0.8996667 0.5004233 0.2989493 1.0935090 0.5941743 ) 
while the true parameter is: 

9* = ( -0.80 0.90 0.50 0.30 1.10 0.60 ) . 

4 Mixture of intervational and observational data 

4.1 Case of a single intervention experiment 

We assume now that we perform an intervention on a subset J C T = {1, ... ,p} of variables by artificially 
setting the level of the corresponding variables to a value: do(Xj = xj). The corresponding model is obtained 
by assuming that all Wij — for e £ and j e J; we denote the corresponding matrix Wj. We also 

assume that the variables Xj for j e J are fully deterministic. The resulting model is hence Gaussian once again: 
X x \do{Xj = xj) ~ N{hj{xj), Sj) with 

V>j{xj) = vj(xj)Lj, Sj = L^diag{a 2 )DjLj 
where Dj = J2j$j e J e j i s a diagonal matrix with at J' positions and 1 elsewere, and with 



{Z 3 else^ and L J = (I-W J )-=I + W J + ... + Wr. 
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4.2 Maximum likelihood estimator 
4.2.1 Likelihood 

We consider N data generated under x k = (x\ 1 . . . , x k ) (1 ^ k ^ N) with intervention on Jk (Jk = means no 
intervention). We denote by K.j = {k,j Jk}, and by Nj = \fCj \ its cardinal. The log-likelihood of the model 
can then be written as: 



£(m, a, w) 



log(27r) 



E^-E N i lo ^) - J E i E (4 - ** We J - ^) 2 - 



2 ^"T °^ fee/c 3 



Proof. This is mainly due to the fact that for any intervention set J we have W jej = WeJ for all j £ J. 
Considering the derivative with respect to rrij we get for all j such that Nj > 0: 



rrij 



k£Kj 



which can be plugged into the likelihood expression to get: 

E * - E - J E i E (tf J - y k '^f 



!(a, w) 



2 ^ crj 

j 3 kelCj 



where for (fc, j) such as fc e Kj we have: 



k L 



k'EJCj 



4.2.2 Estimators 

It can be shown that w may be estimated by solving the following linear system: 

■ k ' j - k ' j forall(i,i) e£. 



E Wi, 'j E 

i',(i',j)e£ keKj 



k,j k,j 



E Vi "yj 



k£Kj 



Note that the system might be degenerate if the intervention design gives no insight on some parameters. 
It is hence finally possible to obtain an estimator of a through: 



^ = ^E(^-^ e j) 2 . 



k^tCj 



4.2.3 Toy example (continued) 

Let us consider the following design: J\ = {1} with x\ = —0.5, Ji = {2} with x\ = 0.5, J3 
x\ = 0.1, Ji = {1,2} with x\ = —1.5 and x\ = 2.5, and J§ = $ (no intervention). 
Here is a sample of size N = 5 drawn from our toy-example model: 

/ -0.50000000 0.9391031 0.7665494 \ 

0.47655556 0.5000000 1.4537910 

x= 0.09892252 1.2963643 0.1000000 

-1.50000000 2.5000000 0.3326028 

\ 0.36614988 1.1787898 1.9014714 / 



{3} with 



after centering we get: 



V = 



( -0.81387599 -0.0526149 -0.3852047 \ 

0.16267957 -0.4917180 0.3020369 

-0.21495346 0.3046462 -1.0517541 

-1.81387599 1.5082820 -0.8191513 

V 0.05227389 0.1870718 0.7497173 / 



/ -0.4883575 
0.4881981 
0.1105651 

-1.4883575 1.36191426 

\ 0.3777924 0.04070409 



0.19898261 -0.1561242 \ 
0.63808574 0.5311174 
-0.8226736 
-0.5900708 
0.9787978 / 



0.15827852 



6 



/ -0.2106764 -0.34037011 -0.3470542 \ 

0.7658792 -0.77947324 0.3401873 

0.3882462 0.01689102 -1.0136037 

-1.2106764 1.22052676 -0.7810009 

\ 0.6554735 -0.10068341 0.7878678 J 



We get w by solving: 



0.3934448 0.000000 0.000000 \ / w h2 
0.0000000 2.526338 -2.068933 u) M 
0.0000000 -2.068933 2.223253 / V w 23 



0.1300524 
1.7956242 
-1.1795977 



which gives (reference value in parentheses): 

u>i, 2 = 0.3305481 {w* 12 = -0.8) u> 1)3 = 1.1612114 (w* 13 = 0.9) w 2 , 3 = 0.5500366 (w^ 3 = 0.5). 
We then finally obtain: 

<7i = 0.15853727 (crj = 0.3) a 2 = 0.08815595 {a* 2 = 1.1) a 3 = 0.08745639 (a* 3 = 0.6) 
and the nuisance parameter: 

mi = 0.3138760 (m* = 0.5) rh 2 = 1.1419342 (m 2 = 0.2) m 3 = 0.7458125 (m* 3 = 0.7). 



4.3 Fisher information 
4.3.1 Hessian of I 

The (non-zero) second order derivatives of I are given by: 



if- 



—(<t,w) = -/-2 E v^ 3 y k ^ 

' i u A , 



d 2 



3 keiCj 



dwijdaj 



(a, w) 



G 



E ^'(yf-y^'wej). 



3 k£Kj 



1 k£Kj 



4.3.2 Distribution of y k ' j 

For all k, let us adopt the following notation: = W j k ,l>k = L , v>k = v j k , fi k = fi j , and = X . 
We can then rewrite Equation (|4) as: 



y fej = x k - 



- E 



Ni - 1 



N, ~ X N. 



- E 



„k' 



k' £Kj ,k'^k 



with x ~ A/" (/x fc ; Sfe), from which we derive that: 



77" 

3 fc'e/Cj 



AT| 



J k'eKj,k'jtk 



\ 



and W>T -A'(0;^-ia-] 



4.3.3 Information 

The Fisher information matrix I(cr, w) can therefore be written as: 

dH 



e m J(a,w)e w _ l 



-E 



dwijdwi' j 



(cr, to) 



i E 



J fce/c,- 



m ■ J m., 



1 / N ] + Nj - 1 



3 k£K 



E s fc + E m 



k , 7 /s , 7 



k£JCj 



7 



M I(<7, w)el 



-E 



d 2 e 



dwijdaj 



(a, w) 



-y 

3 fc6/C, 



E 



.'/"'.'/'•'(J 



W)eJ 



and finally: 



e a -I(a,w)el 



-E 



d 2 l 
dwijduj 



(a, w) 



a* 2^ 



E 



N 



N 



2Ni 



4.3.4 Toy example (continued) 

We consider the same intervention design as before, except that each condition is repeated 40 times. Let us consider 
the following design: J k = {1} with x\ = -0.5 for k = 1 ... 40, J k = {2} with x\ = 0.5 for k = 41 . . . 80, 
Jk = {3} with xl = 0.1 for k = 81 . . . 120, J k = {1, 2} with = -1.5 and x\ = 2.5 for k = 121 . . . 160, and 
Jfc = (no intervention) for k = 161 . . . 200. We thus have /Ci = {41, . . . , 120, 161, . . . , 200} with N x = 120, 
K 2 = {1, • • • , 40, 81, ... , 120, 161, . . . , 200} with N 2 = 120, and K 3 = {1, . . . , 80, 121, . . . , 200} with iV 3 = 
160. 

With this design, we obtain the following Fisher information matrix for (wi t 2, ^1,3 5 $1,2, <5i, 02, 03): 



/ 27.93981 






V 





325.4313 
291.2836 









291.2836 
541.3569 











2633.3333 











195.8678 





\ 





880.5556 / 



which is consistent with the inverse of the empirical covariance matrix (sample size 2000): 



27.17689131 

-1.99150498 
-0.29200154 
-1.06047181 
-0.42604171 
-0.03815717 



-1.991505 
311.458993 
-278.425060 

2.672117 
-6.198404 
4.454986 



-0.2920015 
278.4250602 
519.5078012 

-9.8669528 
7.8182677 
-6.6410163 



-1.060472 
2.672117 
-9.866953 
2708.375488 

-10.567025 
13.185252 



-0.4260417 
-6.1984037 
7.8182677 
-10.5670252 
194.9267470 
-2.8052306 



-0.03815717 \ 

4.45498604 
-6.64101633 
13.18525152 
-2.80523059 
901.29551952 / 



5 Conclusion 

Joint causal network inference from a mixture of observational and intervention transcriptomic data is a very im- 
portant and challenging research question. In this technical note, we provided an explicit formula for the likelihood 
function in the context of Gaussian Bayesian networks under any complex intervention design, as well as its an- 
alytical maximization. For an unknown graph structure with a known parental node order, it is therefore possible 
to directly estimate the causal effects. A crucial next step will be to propose an algorithm to obtain the optimal 
parental order. To this end, we envisage the use of a Mallow's ||8][9) proposal distribution in an empirical Bayesian 
algorithm. 

The choice of optimal experimental intervention designs is an important practical question for biologists plan- 
ning future gene knock-out experiments. Recently, Hauser and Biihlmann (2) proposed two strategies for the 
choice of optimal interventions for Gaussian Bayesian networks. The first is a greedy approach using single-vertex 
interventions that maximize the number of edges that can be oriented after each intervention, and the second yields 
a minimum set of targets of arbitrary size that guarantee full identifiability. Future research will be needed to de- 
termine whether the optimal knock-outs to be performed could alternatively be chosen by evaluating the amount of 
information potentially contributed by each possible intervention via the Fisher information matrix. We note that 
the derivation of the Fisher information matrix is not trivial, especially in the case of a mixture of observational 
and intervention data; in this technical note, we provided formulae for the calculation of the Fisher information, 
providing an opportunity for future research concerning optimal experimental intervention designs. 
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